A/B Testing

Tutorial Source:

https://www.business-science.io/business/2019/03/11/ab-testing-machine-learning.html

A/B Testing

• aka split testing or bucket testing
• in marketing, A/B testing enables us to determine whether changes in landing pages, popup forms, article titles, and other digital marketing decisions improve conversion rates and ultimately customer purchasing favor
• a successful A/B testing strategy can lead to massive gains - more satisfied users, more engagement, and more sales

Method

2 tests are run in parallel: 1. Treatment Group (Group A) - group exposed to the new web page, popup form, etc 2. Control Group (Group B) - group experiences no change from the current setup - goal is to compare the conversion rates of the two groups using statistical inference

ML - A/B Testing

• tidyverse and tidyquant: These are the core data manipulation and visualization packages. We’ll mainly be using dplyr for data manipulation, ggplot2 for data visualization, and tidyquant themes for business reporting.
• parsnip, rsample, recipes, and yardstick: These are the tidyverse modeling packages. The parsnip package is an amazing tool that connects to the main machine learning algorithms. We teach parsnip in-depth (44 lessons, 5 hours of video) in Business Analysis with R, Week 6, Part 2 - Machine Learning (Regression).
• rpart, rpart.plot, and xgboost: These are the modeling libraries that we’ll connect to through the parsnip interface.

# Core packages
library(tidyverse)
library(tidyquant)

# Modeling packages
library(parsnip)
library(recipes)
library(rsample)
library(yardstick)
library(broom)

# Connector packages
library(rpart)
library(rpart.plot)
library(xgboost)

# Other packages
library(data.table)
library(plotly)

control_tbl <- fread("control_data.csv")
experiment_tbl <- fread("experiment_data.csv")

Investigate Data

head(control_tbl) 

glimpse(control_tbl)

Rows: 37
Columns: 5
$ Date        <chr> "Sat, Oct 11", "Sun, Oct 12", "Mon, Oct 13", "Tue, Oct 14", "Wed, Oct 15", "Thu, Oct 16", "Fri, Oct…
$ Pageviews   <int> 7723, 9102, 10511, 9871, 10014, 9670, 9008, 7434, 8459, 10667, 10660, 9947, 8324, 9434, 8687, 8896,…
$ Clicks      <int> 687, 779, 909, 836, 837, 823, 748, 632, 691, 861, 867, 838, 665, 673, 691, 708, 759, 736, 739, 734,…
$ Enrollments <int> 134, 147, 167, 156, 163, 138, 146, 110, 131, 165, 196, 162, 127, 220, 176, 161, 233, 154, 196, 167,…
$ Payments    <int> 70, 70, 95, 105, 64, 82, 76, 70, 60, 97, 105, 92, 56, 122, 128, 104, 124, 91, 86, 75, 101, 93, 67, …

glimpse(experiment_tbl)

Rows: 37
Columns: 5
$ Date        <chr> "Sat, Oct 11", "Sun, Oct 12", "Mon, Oct 13", "Tue, Oct 14", "Wed, Oct 15", "Thu, Oct 16", "Fri, Oct…
$ Pageviews   <int> 7716, 9288, 10480, 9867, 9793, 9500, 9088, 7664, 8434, 10496, 10551, 9737, 8176, 9402, 8669, 8881, …
$ Clicks      <int> 686, 785, 884, 827, 832, 788, 780, 652, 697, 860, 864, 801, 642, 697, 669, 693, 771, 736, 727, 728,…
$ Enrollments <int> 105, 116, 145, 138, 140, 129, 127, 94, 120, 153, 143, 128, 122, 194, 127, 153, 213, 162, 201, 207, …
$ Payments    <int> 34, 91, 79, 92, 94, 61, 44, 62, 77, 98, 71, 70, 68, 94, 81, 101, 119, 120, 96, 67, 123, 100, 103, N…

Conclusions: - data is in character format; need to convert to date - payment is an outcome of enrollments so this should be removed

Check for Missing Data

control_tbl %>% 
  map_df(~ sum(is.na(.))) %>% 
  gather(key = "feature", value = "missing_count") %>% 
  arrange(desc(missing_count))

experiment_tbl %>% 
  map_df(~ sum(is.na(.))) %>% 
  gather(key = "feature", value = "missing_count") %>% 
  arrange(desc(missing_count))

control_tbl %>% 
  filter(is.na(Enrollments))

• we don’t have Enrollment information from November 3rd on. We will need to remove these observations

Format Data

• Combine the control_tbl and experiment_tbl, adding an “id” column indicating if the data was part of the experiment or not
• Add a “row_id” column to help for tracking which rows are selected for training and testing in the modeling section
• Create a “Day of Week” feature from the “Date” column
• Drop the unnecessary “Date” column and the “Payments” column
• Handle the missing data (NA) by removing these rows.
• Shuffle the rows to mix the data up for learning
• Reorganize the columns

set.seed(123)

data_formatted_tbl <- control_tbl %>% 
  
  # combine with experiment data
  bind_rows(experiment_tbl, .id = "Experiment") %>% 
  mutate(Experiment = as.numeric(Experiment) - 1) %>% 
  
  # add row id
  mutate(row_id = row_number()) %>% 
  
  # create a day of week feature
  mutate(DOW = str_sub(Date, start = 1, end = 3) %>% 
           factor(levels = c("Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"))
         ) %>% 
  select(-Date, -Payments) %>% 
  
  # remove missing data
  filter(!is.na(Enrollments)) %>% 
  
  # shuffle the data (note that set.seed is used to make reproducible)
  slice_sample(prop = 1) %>% 
  
  # reorganize columns
  select(row_id, Enrollments, Experiment, everything())

glimpse(data_formatted_tbl)

Rows: 46
Columns: 6
$ row_id      <int> 45, 15, 14, 3, 56, 51, 58, 39, 40, 41, 5, 55, 42, 9, 43, 57, 8, 52, 7, 10, 47, 19, 4, 54, 17, 11, 4…
$ Enrollments <int> 94, 176, 220, 167, 201, 194, 182, 116, 145, 138, 163, 162, 140, 131, 129, 207, 110, 127, 146, 165, …
$ Experiment  <dbl> 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, …
$ Pageviews   <int> 7664, 8687, 9434, 10511, 9262, 9402, 8715, 9288, 10480, 9867, 10014, 9396, 9793, 8459, 9500, 9308, …
$ Clicks      <int> 652, 691, 673, 909, 727, 697, 722, 785, 884, 827, 837, 736, 832, 691, 788, 728, 632, 669, 748, 861,…
$ DOW         <fct> Sat, Sat, Fri, Mon, Wed, Fri, Fri, Sun, Mon, Tue, Wed, Tue, Wed, Sun, Thu, Thu, Sat, Sat, Fri, Mon,…

Split Data: Training and Testing

set.seed(123)

split_obj <- data_formatted_tbl %>% data.table() %>% 
  initial_split(prop = 0.8, strata = "Experiment")

train_tbl <- training(split_obj)
test_tbl <- testing(split_obj)

glimpse(train_tbl)

Rows: 36
Columns: 6
$ row_id      <int> 15, 3, 5, 9, 8, 19, 4, 17, 11, 13, 20, 23, 18, 21, 6, 1, 2, 16, 45, 56, 51, 58, 39, 41, 55, 42, 43,…
$ Enrollments <int> 176, 167, 163, 131, 110, 196, 156, 233, 196, 127, 167, 206, 154, 174, 138, 134, 147, 161, 94, 201, …
$ Experiment  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ Pageviews   <int> 8687, 10511, 10014, 8459, 7434, 9327, 9871, 9535, 10660, 8324, 9345, 8836, 9363, 8890, 9670, 7723, …
$ Clicks      <int> 691, 909, 837, 691, 632, 739, 836, 759, 867, 665, 734, 693, 736, 706, 823, 687, 779, 708, 652, 727,…
$ DOW         <fct> Sat, Mon, Wed, Sun, Sat, Wed, Tue, Mon, Tue, Thu, Thu, Sun, Tue, Fri, Thu, Sat, Sun, Sun, Sat, Wed,…

glimpse(test_tbl)

Rows: 10
Columns: 6
$ row_id      <int> 14, 40, 52, 7, 10, 12, 59, 49, 22, 60
$ Enrollments <int> 220, 145, 127, 146, 165, 162, 142, 128, 156, 182
$ Experiment  <dbl> 0, 1, 1, 0, 0, 0, 1, 1, 0, 1
$ Pageviews   <int> 9434, 10480, 8669, 9008, 10667, 9947, 8448, 9737, 8460, 8836
$ Clicks      <int> 673, 884, 669, 748, 861, 838, 695, 801, 681, 724
$ DOW         <fct> Fri, Mon, Sat, Fri, Mon, Wed, Sat, Wed, Sat, Sun

Implement ML Algorithms

• implement 3 modeling approaches:

1. Linear Regression - Linear, Explainable (Baseline)
2. Decision Tree - Pros: Non-Linear, Explainable; Cons: Lower Performance
3. XGBoost - Pros: Non-Linear, High Performance; Cons: Less Explainable

Linear Regression (Baseline)

• create a model using the training data and predict on the test data

model_01_lm <- linear_reg("regression") %>% 
  set_engine("lm") %>% 
  fit(Enrollments ~ ., data = train_tbl %>% select(-row_id))

pred_01_lm <- model_01_lm %>%
  predict(new_data = test_tbl) %>%
  bind_cols(test_tbl %>% select(Enrollments)) 

pred_01_lm %>% 
  metrics(truth = Enrollments, estimate = .pred) %>% 
  knitr::kable()

.metric	.estimator	.estimate
rmse	standard	32.1813991
rsq	standard	0.3444336
mae	standard	27.7148176

pred_01_lm %>% 
  mutate(observation = row_number() %>% as.character()) %>% 
  gather(key = "key", value = "value", -observation, factor_key = TRUE) %>% 
  plot_ly(x = ~observation, y = ~value, color = ~key, type = "scatter", mode = "markers") %>%
  layout(title = "Prediction vs Annual Enrollments: Model 1 - Linear Regression")

What’s driving this model?

Use the tidy() function from the broom package to help. This gets us the model estimates. We can arrange by “p.value” to get an idea of how important the model terms are. Clicks, Pageviews, and Experiment are judged strong predictors with a p-value less than 0.05. However, we want to try out other modeling techniques to judge this. We note that the coefficient of Experiment is -17.6, and because the term is binary (0 or 1) this can be interpreted as decreasing Enrollments by -17.6 per day when the Experiment is run.

linear_regression_model_terms_tbl <- model_01_lm$fit %>% 
  tidy() %>% 
  arrange(p.value) %>% 
  mutate(term = as.factor(term) %>% fct_rev())

linear_regression_model_terms_tbl %>% knitr::kable()

term	estimate	std.error	statistic	p.value
Clicks	-0.6709381	0.1351726	-4.9635647	0.0000370
Pageviews	0.0692229	0.0154320	4.4856713	0.0001306
Experiment	-16.8809060	7.3964356	-2.2823029	0.0308955
DOWMon	31.1476094	17.1397658	1.8172716	0.0807195
DOWFri	18.0185100	14.0265834	1.2845972	0.2102610
DOWWed	18.3067550	15.1370069	1.2094039	0.2373904
DOWSat	10.2618805	16.4042902	0.6255608	0.5370559
DOWThu	-6.6848206	12.3334080	-0.5420092	0.5924275
(Intercept)	26.5108797	77.5271192	0.3419562	0.7351325
DOWTue	-4.6330479	15.4975121	-0.2989543	0.7673515

Visualizing the features show that clicks, pageviews and experiment are the most important features and are likely to be the best predictors

linear_regression_model_terms_tbl %>% 
  arrange(desc(p.value), term) %>% 
  plot_ly(x = ~p.value, y = ~term, type = "scatter", mode = "markers") %>% 
  add_lines(x = ~0.05) %>% 
  layout(title = "Feature Importance: Model 1 - Linear Regression")

Set-up Helper Functions

calc_metrics <- function(model, new_data) {
  model %>%
    predict(new_data = new_data) %>% 
    bind_cols(new_data %>% select(Enrollments)) %>% 
    metrics(truth = Enrollments, 
            estimate = .pred)
}

plot_predictions <- function(model, new_data) {
  predict(model, new_data) %>% 
    bind_cols(new_data %>% select(Enrollments)) %>% 
    mutate(observation = row_number() %>% as.character()) %>% 
    gather(key = "key", value = "value", -observation, factor_key = TRUE) %>% 
    plot_ly(x = ~observation, y = ~value, color = ~key, type = "scatter", mode = "markers")
}

Decision Trees

Decision Trees are excellent models that can pick up on non-linearities and often make very informative models that compliment linear models by providing a different way of viewing the problem. We can implement a decision tree with decision_tree(). We’ll set the engine to “rpart”, a popular decision tree package. There are a few key tunable parameters: - cost_complexity: A cutoff for model splitting based on increase in explainability - tree_depth: The max tree depth - min_n: The minimum number of observations in terminal (leaf) nodes The parameters selected for the model were determined using 5-fold cross validation to prevent over-fitting.

model_02_decision_tree <- decision_tree(
    mode = "regression", 
    cost_complexity = 0.001, 
    tree_depth = 5, 
    min_n = 4) %>% 
  set_engine("rpart") %>% 
  fit(Enrollments ~ ., data = train_tbl %>% select(-row_id))

Next, calculate the metrics on this model using our helper function, calc_metrics(). The MAE of the predictions is approximately the same as the linear model at +/-19 Enrollments per day.

model_02_decision_tree %>% 
  calc_metrics(test_tbl) %>% 
  knitr::kable()

.metric	.estimator	.estimate
rmse	standard	22.71472
rsq	standard	0.38662
mae	standard	18.08333

model_02_decision_tree %>% 
  plot_predictions(test_tbl) %>% 
  layout(title = "Prediction vs Actual Enrollments: Decision Tree")

Warning in RColorBrewer::brewer.pal(N, "Set2") :
  minimal value for n is 3, returning requested palette with 3 different levels

Warning in RColorBrewer::brewer.pal(N, "Set2") :
  minimal value for n is 3, returning requested palette with 3 different levels

Warning in RColorBrewer::brewer.pal(N, "Set2") :
  minimal value for n is 3, returning requested palette with 3 different levels

Warning in RColorBrewer::brewer.pal(N, "Set2") :
  minimal value for n is 3, returning requested palette with 3 different levels

Finally, use rpart.plot() to visualize the decision tree rules. Note that we need to extract the underlying “rpart” model from the parsnip model object using the model_02_decision_tree$fit

model_02_decision_tree$fit %>% 
  rpart.plot(roundint = FALSE, 
             cex = 0.8, 
             fallen.leaves = TRUE, 
             extra = 101, 
             main = "Model 02: Decision Tree")

Interpreting the decision tree is straightforward: Each decision is a rule, and Yes is to the left, No is to the right. The top features are the most important to the model (“Pageviews” and “Clicks”). The decision tree shows that “Experiment” is involved in the decision rules. The rules indicate a when Experiment >= 0.5, there is a drop in enrollments.

Key Points:

Our new model has roughly the same accuracy to +/-19 enrollments (MAE) as the linear regression model.

Experiment shows up towards the bottom of the tree. The rules indicate a when Experiment >= 0.5, there is a drop in enrollments.

XGBoost

Several key tuning parameters include: - mtry: The number of predictors that will be randomly sampled at each split when creating the tree models. - trees: The number of trees contained in the ensemble. - min_n: The minimum number of data points in a node that are required for the node to be split further. - tree_depth: The maximum depth of the tree (i.e. number of splits). - learn_rate: The rate at which the boosting algorithm adapts from iteration-to-iteration. - loss_reduction: The reduction in the loss function required to split further. - sample_size: The amount of data exposed to the fitting routine.

The parameters selected for the model were determined using 5-fold cross validation to prevent over-fitting.

set.seed(123)

model_03_xgboost <- boost_tree(
    mode = "regression", 
    mtry = 100, 
    trees = 1000, 
    min_n = 8, 
    tree_depth = 6, 
    learn_rate = 0.2, 
    loss_reduction = 0.01, 
    sample_size = 1) %>% 
  set_engine("xgboost") %>% 
  fit(Enrollments ~ ., data = train_tbl %>% select(-row_id))

model_03_xgboost %>% 
  calc_metrics(test_tbl) %>% 
  knitr::kable()

.metric	.estimator	.estimate
rmse	standard	20.5911658
rsq	standard	0.4683664
mae	standard	17.5331612

model_03_xgboost %>% 
  plot_predictions(test_tbl) %>% 
  layout(title = "Prediction vs Actual Enrollments: XGBoost")

Warning in RColorBrewer::brewer.pal(N, "Set2") :
  minimal value for n is 3, returning requested palette with 3 different levels

Warning in RColorBrewer::brewer.pal(N, "Set2") :
  minimal value for n is 3, returning requested palette with 3 different levels

Warning in RColorBrewer::brewer.pal(N, "Set2") :
  minimal value for n is 3, returning requested palette with 3 different levels

Warning in RColorBrewer::brewer.pal(N, "Set2") :
  minimal value for n is 3, returning requested palette with 3 different levels

Feature Importance

xgboost_feature_importance_tbl <- model_03_xgboost$fit %>% 
  xgb.importance(model = .) %>% 
  as_tibble() %>% 
  mutate(Feature = as_factor(Feature) %>% fct_rev())

xgboost_feature_importance_tbl %>% knitr::kable()

Feature	Gain	Cover	Frequency
Pageviews	0.5551216	0.5979430	0.5762712
Clicks	0.3824186	0.2145172	0.2364613
Experiment	0.0624599	0.1875398	0.1872675

xgboost_feature_importance_tbl %>% 
  mutate(Label = paste0(round(Gain*100, 1)), "%") %>% 
  plot_ly(x = ~Gain, y = ~Feature, type = "scatter", mode = "markers", name = ~Label) %>% 
  layout(title = "XGBoost Feature Importance")

The information gain is 93% from Pageviews and Clicks combined. Experiment has about a 7% contribution to information gain, indicating it’s still predictive (just not nearly as much as Pageviews). This tells a story that if Enrollments are critical, Udacity should focus on getting Pageviews.

Key Points: - The XGBoost model error has dropped to +/-11 Enrollments. - The XGBoost shows that Experiment provides an information gain of 7% - The XGBoost model tells a story that Udacity should be focusing on Page Views and secondarily Clicks to maintain or increase Enrollments. The features drive the system.

Business Conclusions

There are several key benefits to performing A/B Testing using Machine Learning. These include:

Understanding the Complex System - We discovered that the system is driven by Pageviews and Clicks. Statistical Inference would not have identified these drivers. Machine Learning did.

Providing a direction and magnitude of the experiment - We saw that Experiment = 1 drops enrollments by -17.6 Enrollments Per Day in the Linear Regression. We saw similar drops in the Decision Tree rules. Statistical inference would not have identified magnitude and direction. Only whether or not the Experiment had an effect.

What Should Udacity Do?

If Udacity wants to maximimize enrollments, it should focus on increasing Page Views from qualified candidates. Page Views is the most important feature in 2 of 3 models.

If Udacity wants alert people of the time commitment, the additional popup form is expected to decrease the number of enrollments. The negative impact can be seen in the decision tree (when Experiment <= 0.5, Enrollments go down) and in the linear regression model term (-17.6 Enrollments when Experiment = 1). Is this OK? It depends on what Udacity’s goals are.

Cross Validation and Improving Modeling Performance

Two important further considerations when implementing an A/B Test using Machine Learning are: 1. How to Improve Modeling Performance 2. The need for Cross-Validation for Tuning Model Parameters

Improving Modeling Performance

• run analysis on unaggregated data (data in this exercise was aggregated)
• run analysis on individual customer data to determine probability on an individual customer enrolling
• include good features; customer-related features not included in this data set

Cross-Validation for Tuning Models

• In practice, we need to perform cross-validation to prevent the models from being tuned to the test data set.

Using caret